Run UMAP on CD4+ events which express a COMPASS subset. Repeat for CD8 events.
The question being asked is “What are the memory and activation profiles of Ag-specific T cells?”

This time around, don’t sample the events.
Don’t stratify by groups, but rather color the sub-localization of the different markers.
Color by Cohort and Antigen (S1, S2, NCAP, VEMP, including DMSO)
Also color by
- Degree of functionality
- Cytokine - CD45RA
- CCR7
- HLA-DR
- CD38
Boxplots?

library(openCyto)
library(CytoML)
library(flowCore)
library(flowWorkspace)
library(here)
library(tidyverse)
library(uwot)
library(ggplot2)
library(scales)
library(patchwork)
library(hues)
library(RColorBrewer)
library(ggrepel)
library(ggpubr)
library(tidyselect)
library(ggrastr)
source(here::here("scripts/20200604_Helper_Functions.R")) # for distributeEvents() and sampleGatingHierarchy()
date <- 20200815
save_output <- FALSE
rerun_dimred <- FALSE

gs <- load_gs(here::here("out/GatingSets/20200815_HAARVI_ICS_GatingSet_AllBatches_with_COMPASS_Subsets_R4.0.3"))

gs2 <- subset(gs, !(`SAMPLE ID` %in% c("37C", "BWT23", "116C", "BWT22")) &
                !(`SAMPLE ID` == "551432" & STIM == "Spike 2"))
dput(gh_get_pop_paths(gs2))
## c("root", "/Time", "/Time/LD-3+", "/Time/LD-3+/1419-3+", "/Time/LD-3+/1419-3+/S", 
## "/Time/LD-3+/1419-3+/S/Lymph", "/Time/LD-3+/1419-3+/S/Lymph/4+", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/107a", "/Time/LD-3+/1419-3+/S/Lymph/4+/154", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CCR7+", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD45RA+", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL2", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL4513", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/TNF", "/Time/LD-3+/1419-3+/S/Lymph/8+", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/CD38+", 
## "/Time/LD-3+/1419-3+/S/Lymph/HLADR+", "/Time/LD-3+/1419-3+/S/Lymph/NOT4+", 
## "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/107a", "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/154", 
## "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/CCR7+", "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/CD45RA+", 
## "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/IL2", 
## "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/IL4513", 
## "/Time/LD-3+/1419-3+/S/Lymph/NOT4+/TNF", "/Time/LD-3+/1419-3+/S/Lymph/8+/107a", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/154", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL2", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL4513", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/TNF", "/Time/LD-3+/1419-3+/S/Lymph/8+/CCR7+", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD45RA+", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_IL2_AND_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_IL17_AND_NOT_IL2_AND_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_NOT_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_NOT_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_NOT_IFNG_AND_NOT_IL17_AND_IL2_AND_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_NOT_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_107a_AND_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets", "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_NOT_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_107a_AND_NOT_154_AND_NOT_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_NOT_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_NOT_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_107a_AND_NOT_154_AND_IFNG_AND_NOT_IL17_AND_NOT_IL2_AND_NOT_IL4513_AND_TNF", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets/Naive", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets/TCM", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets/TEMRA", 
## "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets/TEM", "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets/Naive", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets/TCM", "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets/TEMRA", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets/TEM", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets/HLADR+CD38+", 
## "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets/HLADR+CD38+"
## )

CD4

Extract mfi data

cd4_gates_for_dimred <- c(
  "/Time/LD-3+/1419-3+/S/Lymph/4+/107a", "/Time/LD-3+/1419-3+/S/Lymph/4+/154", 
  "/Time/LD-3+/1419-3+/S/Lymph/4+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL2", 
  "/Time/LD-3+/1419-3+/S/Lymph/4+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL4513", 
  "/Time/LD-3+/1419-3+/S/Lymph/4+/TNF",
  "/Time/LD-3+/1419-3+/S/Lymph/4+/CCR7+", "/Time/LD-3+/1419-3+/S/Lymph/4+/CD45RA+",
  "/Time/LD-3+/1419-3+/S/Lymph/CD38+", "/Time/LD-3+/1419-3+/S/Lymph/HLADR+")
cd4_cytokine_gates <- c("/Time/LD-3+/1419-3+/S/Lymph/4+/107a", "/Time/LD-3+/1419-3+/S/Lymph/4+/154", 
                        "/Time/LD-3+/1419-3+/S/Lymph/4+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL2", 
                        "/Time/LD-3+/1419-3+/S/Lymph/4+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/4+/IL4513", 
                        "/Time/LD-3+/1419-3+/S/Lymph/4+/TNF")
cd4_compass_subsets_parentGate <- "/Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets"

pop_counts <- pData(gs2) %>% 
  left_join(gs_pop_get_count_fast(gs2, subpopulations = cd4_compass_subsets_parentGate),
            by = c("rowname" = "name")) %>% 
  dplyr::rename(CD4_COMPASS_Subsets = Count) %>% 
  dplyr::select(rowname, Batch, "SAMPLE ID", STIM, Cohort, CD4_COMPASS_Subsets) %>% 
  dplyr::filter(!(Cohort %in% c(NA, "Healthy control", "Healthy control 2017-2018")) & STIM != "SEB")

Keep in mind that there is lopsided patient and group representation simply due to not sampling:

cd4_compass_subsets_sampleSizes_4plot <- pop_counts %>%
  mutate(Cohort = factor(Cohort,
                         levels = c("Non-hospitalized", "Hospitalized"),
                         labels = c("Conv\nNon-Hosp", "Conv\nHosp")))
ggplot(cd4_compass_subsets_sampleSizes_4plot,
       aes(factor(Cohort), CD4_COMPASS_Subsets)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.15, height = 0) +
  theme_bw(base_size=20) +
  labs(title="ICS CD4 UMAP patient representation",
       y="CD4+ COMPASS Subset+ Events\n for Dimensionality Reduction\n(not sampled)") +
  facet_grid(Batch ~ STIM) +
  theme(axis.title.x = element_blank())

# Extract data for dimensionality reduction (not actually sampling)
call_sampleGatingHierarchy_for_cd4 <- function(currentSampleName) {
  # print(sprintf("Sampling data from %s", currentSampleName))
  sampleGatingHierarchy(gs2[[currentSampleName]], cd4_compass_subsets_parentGate, n = NULL, otherGates = cd4_gates_for_dimred)
}
cd4_compass_subsets_data <- map_dfr(pop_counts$rowname, call_sampleGatingHierarchy_for_cd4)
dim(cd4_compass_subsets_data)
## [1] 123642     54
knitr::kable(head(cd4_compass_subsets_data))
Days symptom onset to visit 1 Sex WELL ID Pair ID Hispanic? STIM name Cell count Race Cohort Age Race_v2 rowname $DATE Batch Sample ID SAMPLE ID PLATE NAME Collection date EXPERIMENT NAME filename /Time/LD-3+/1419-3+/S/Lymph/4+/CD4_COMPASS_Subsets /Time/LD-3+/1419-3+/S/Lymph/4+/107a /Time/LD-3+/1419-3+/S/Lymph/4+/154 /Time/LD-3+/1419-3+/S/Lymph/4+/IFNG /Time/LD-3+/1419-3+/S/Lymph/4+/IL2 /Time/LD-3+/1419-3+/S/Lymph/4+/IL17 /Time/LD-3+/1419-3+/S/Lymph/4+/IL4513 /Time/LD-3+/1419-3+/S/Lymph/4+/TNF /Time/LD-3+/1419-3+/S/Lymph/4+/CCR7+ /Time/LD-3+/1419-3+/S/Lymph/4+/CD45RA+ /Time/LD-3+/1419-3+/S/Lymph/CD38+ /Time/LD-3+/1419-3+/S/Lymph/HLADR+ Time FSC-A FSC-H SSC-A SSC-H CD8b TNFa CD107a CD154 CD3 ECD IL2 CD4 IL17a IL4/5/13 CD14/CD19 CCR7 CD38 L/D IFNg CD45RA HLADR
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 1 0 1 0 0 0 15.439 122568.24 100216 31674.96 30776 1027.0212 758.3488 158.9819 383.4604 3007.214 993.8220 1901.014 988.5380 1611.0039 663.5084 1872.624 1194.567 1324.692 432.7889 1398.416 967.1333
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 1 0 1 0 0 1 1 0 0 0 15.485 142859.48 123614 32029.05 29130 286.6093 2734.2783 1114.5400 2350.0759 2810.713 3142.3025 1715.415 1170.4102 833.1918 492.7284 1535.650 1198.652 1217.047 1484.1797 1539.542 173.3101
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 0 1 0 0 0 0 15.562 120323.20 105313 36177.21 35201 645.6627 1941.5667 786.9394 1029.9525 2903.400 1293.5875 1862.381 855.8306 635.2674 954.0318 1328.853 1188.295 1422.533 939.9696 1395.431 699.8994
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 0 1 1 0 0 0 16.042 65290.08 48154 27090.93 25210 349.0576 1358.4554 504.7423 1294.9139 2929.793 734.6826 2022.920 1008.4814 305.6316 765.4292 1983.099 1334.761 1229.228 962.7062 1490.717 853.8486
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 0 1 0 0 0 0 16.486 66916.48 55401 17046.78 16637 664.9983 1958.0522 787.3279 1061.2051 2986.041 1464.6097 1927.694 1154.1066 385.6267 740.7294 1373.256 1464.787 1084.778 994.0617 1497.823 946.3115
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 1 0 1 1 0 0 16.574 76877.80 62924 16660.50 16012 163.2860 760.8881 518.0920 1345.3741 3152.756 393.7231 1930.848 1694.5745 2547.2505 812.1132 2584.446 1393.654 1234.604 537.8035 2242.353 229.1811

Is there maybe one cytokine that is dominating the entire sample?
If one cytokine has very high background expression (and a generous gate), it could be gated positive in a lot of events.
The high number of events expressing this cytokine could lead to it dominating the data, so that most sampled events are positive for this noisy cytokine. It would drown out real signal from other cytokines.

cytokine_dominance <- cd4_compass_subsets_data %>% 
  group_by(Batch) %>% 
  summarise_at(cd4_cytokine_gates, sum) %>%
  t() %>%
  as.data.frame() %>%
  set_names(c("B1", "B2", "B3")) %>%
  rownames_to_column("Cytokine_Gate") %>%
  dplyr::filter(Cytokine_Gate != "Batch")
knitr::kable(cytokine_dominance)
Cytokine_Gate B1 B2 B3
/Time/LD-3+/1419-3+/S/Lymph/4+/107a 10793 16208 13801
/Time/LD-3+/1419-3+/S/Lymph/4+/154 10689 12468 12604
/Time/LD-3+/1419-3+/S/Lymph/4+/IFNG 7405 4551 5493
/Time/LD-3+/1419-3+/S/Lymph/4+/IL2 10368 8386 9887
/Time/LD-3+/1419-3+/S/Lymph/4+/IL17 618 781 1392
/Time/LD-3+/1419-3+/S/Lymph/4+/IL4513 6537 9705 11946
/Time/LD-3+/1419-3+/S/Lymph/4+/TNF 11258 11151 14956
cytokine_dominance %>%
  pivot_longer(cols = starts_with("B"), names_to = "Batch", values_to = "Events_in_Gate") %>% 
  mutate(Cytokine = sub(".*4\\+\\/(.*)", "\\1", Cytokine_Gate)) %>% 
  ggplot(aes(Cytokine, Events_in_Gate, fill = Cytokine)) +
  theme_bw(base_size=18) +
  geom_bar(stat="identity") +
  facet_grid(. ~ Batch) +
  labs(title = "CD4 Run Cytokine Dominance by Batch")

Perform Dimensionality Reduction

cols_4_dimred <- c("CD3 ECD", "CD8b", "CD4",
                  "TNFa", "CD107a", 
                  "CD154", "IL2", "IL17a", 
                  "IL4/5/13", "IFNg", 
                  "CCR7", "CD45RA",
                  "CD38", "HLADR")
cd4.scaled_dimred_input <- cd4_compass_subsets_data %>% 
  dplyr::select(Batch, all_of(cols_4_dimred)) %>% 
  group_by(Batch) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(data = lapply(data, function(df) {as.data.frame(scale(as.matrix(df)))})) %>% 
  unnest(cols = c(data)) %>% 
  rename_at(vars(all_of(cols_4_dimred)),function(x) paste0(x,".scaled")) %>% 
  dplyr::select(-Batch)

cd4_compass_subsets_data <- cbind(cd4_compass_subsets_data, cd4.scaled_dimred_input)

# UMAP can take a long time, so there is a rerun_dimred switch
if(rerun_dimred) {
  print("Running UMAP")
  set.seed(date)
  print(Sys.time())
  cd4_compass_subsets_dimred_out <- cd4_compass_subsets_data %>% 
    # Run CD3, co-receptor, cytokine, memory, and activation markers through UMAP
    dplyr::select(all_of(paste0(cols_4_dimred, ".scaled"))) %>% 
    uwot::umap(spread = 9, min_dist = 0.02, n_threads = 7)
  print(Sys.time())
  cd4_compass_subsets_w_umap <- cbind(as.data.frame(cd4_compass_subsets_dimred_out) %>% 
    dplyr::rename(x.umap = V1, y.umap = V2),
    cd4_compass_subsets_data)
  if(save_output) {
    saveRDS(cd4_compass_subsets_w_umap, here::here(sprintf("out/UMAP/%s_ICS_CD4_COMPASS_Subsets_UMAP_Unsampled.rds", date)))
  }
} else {
  # Assuming UMAP results are already saved
  print("Loading saved UMAP run")
  cd4_compass_subsets_w_umap <- readRDS(here::here(sprintf("out/UMAP/%s_ICS_CD4_COMPASS_Subsets_UMAP_Unsampled.rds", date)))
}
## [1] "Loading saved UMAP run"

Plot UMAP results

Shuffle data frame rows so e.g. Batch 3 doesn’t dominate foreground

set.seed(date)
cd4_compass_subsets_w_umap <- cd4_compass_subsets_w_umap[sample(nrow(cd4_compass_subsets_w_umap), nrow(cd4_compass_subsets_w_umap)),]
# Arial font setup. Downloaded afms from https://github.com/microsoft/microsoft-r-open/tree/ec3fd89e5fb5794bd8149905c134ad801bb61800
Arial <- Type1Font(family = "Arial",
                   metrics = c(here::here("data/Arial_afm/ArialMT.afm"),
                               here::here("data/Arial_afm/ArialMT-Bold.afm"), 
                               here::here("data/Arial_afm/ArialMT-Italic.afm"),
                               here::here("data/Arial_afm/ArialMT-BoldItalic.afm")))
pdfFonts(Arial = Arial)

boolColorScheme <- c("FALSE" = "#E2E2E2", "TRUE" = "#023FA5")

cd4_compass_subsets_w_umap <- cd4_compass_subsets_w_umap %>% 
  mutate(cytokine_degree = rowSums(dplyr::select(., all_of(cd4_cytokine_gates))))
cd4_compass_subsets_w_umap <- cd4_compass_subsets_w_umap %>% 
  mutate(Cohort = factor(Cohort, levels = c("Non-hospitalized", "Hospitalized")),
         STIM = factor(STIM, levels = c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")))

stim_labs <- c("DMSO", "S1", "S2", "NCAP", "VEMP") 
names(stim_labs) <- c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")
cohort_labs <- c("Conv\nNon-Hosp", "Conv\nHosp")
names(cohort_labs) <- c("Non-hospitalized", "Hospitalized")

base_dimred_plot <- function(currentColumn, pointSize = 0.02, colorScheme = NA) {
  p <- ggplot(cd4_compass_subsets_w_umap, aes(x=x.umap, y=y.umap,
                                       colour=if(currentColumn %in% c("Batch", "SAMPLE ID")) {
                                         factor(!!as.name(currentColumn))
                                         } else {
                                           as.logical(!!as.name(currentColumn))
                                           })) +
    geom_point_rast(shape=20, alpha=0.8, size=pointSize) +
    facet_grid(Cohort ~ STIM, switch="y",
               labeller = labeller(Cohort = cohort_labs, STIM = stim_labs)) +
    theme_bw() +
    theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
          axis.text=element_blank(),
          axis.line=element_blank(),
          axis.ticks=element_blank(),
          axis.title=element_blank(),
          strip.text=element_text(face="bold", size=22),
          panel.grid.major = element_blank(),
          legend.title=element_text(face="bold", size=14),
          strip.text.x = element_text(margin = margin(0.15,0,0.15,0, "cm")))
  if(!anyNA(colorScheme)) {
    p <- p + scale_color_manual(values = colorScheme)
  }
  if(currentColumn %in% c("Batch", "SAMPLE ID")) {
    p <- p + labs(color = currentColumn) +
      guides(colour = guide_legend(override.aes = list(size=10))) +
      theme(legend.title = element_text(size=15),
            legend.text = element_text(size=15),
            legend.position = "bottom") +
      scale_colour_manual(values=as.character(iwanthue(length(unique(cd4_compass_subsets_w_umap[,currentColumn])))))
  } else {
    p <- p + theme(legend.position = "none")
  }
  p
}

Look for Batch and Patient sub-clustering

base_dimred_plot("Batch")

base_dimred_plot("SAMPLE ID")

Cytokine, Memory, and Activation expression localization

Now visualize the cytokine, memory, and activation expression localization

for(cg in cd4_gates_for_dimred) {
  print(base_dimred_plot(cg, colorScheme = boolColorScheme) +
    labs(title = sprintf("CD4+ COMPASS Subset+ UMAP\nColored by %s", sub(".*\\/([^(\\/)]+)", "\\1", cg))))
}

Color by degree

table(cd4_compass_subsets_w_umap$cytokine_degree)
## 
##     1     2     3     4     5 
## 88621 11421 15075  8316   209
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
sc <- scale_colour_gradientn(colours = myPalette(100), limits=c(1, 6))

currentColumn <- "cytokine_degree"
pointSize <- 0.02
ggplot(cd4_compass_subsets_w_umap, aes(x=x.umap, y=y.umap, colour=!!as.name(currentColumn))) +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize) +
  facet_grid(Cohort ~ STIM, switch="y",
             labeller = labeller(Cohort = cohort_labs, STIM = stim_labs)) +
  labs(colour="Cytokine\nDegree",
       title="CD4+ COMPASS Subset+ UMAP\nColored by Cytokine Polyfunctionality") +
  theme_bw() +
  theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        strip.text=element_text(face="bold", size=22),
        panel.grid.major = element_blank(),
        legend.title=element_text(face="bold", size=15),
        legend.text = element_text(size=13),
        strip.text.x = element_text(margin = margin(0.15,0,0.15,0, "cm")),
        legend.position = "bottom") +
  sc

Unstratified UMAP plots

pointSize <- 0.02
# Theme settings
# Expand axis limits so contour lines don't get cut off near borders
gg_themes <- ggplot(cd4_compass_subsets_w_umap, aes(x=x.umap, y=y.umap)) +
  scale_x_continuous(limits=range(cd4_compass_subsets_w_umap$x.umap)*1.1) +
  scale_y_continuous(limits=range(cd4_compass_subsets_w_umap$y.umap)*1.1) +
  theme_bw() +
  theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = margin(1, 0, 0, 0, unit = "pt"))
# Hospitalization status contour plots
hosp_contour_plot <- gg_themes +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
  geom_density_2d(data=cd4_compass_subsets_w_umap %>% dplyr::filter(Cohort == "Hospitalized"),
                  colour="#363636", bins=12) + # "red"
  ggtitle("Hosp")

nonhosp_contour_plot <- gg_themes +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
  geom_density_2d(data=cd4_compass_subsets_w_umap %>% dplyr::filter(Cohort == "Non-hospitalized"),
                  colour="#363636", bins=12) + # "blue"
  ggtitle("Not-Hosp")
# Scaled mfi expression plots for memory, activation, and cytokine markers 
mfi_col_text_4plot <- c("TNFa" = "TNF", "CD107a" = "CD107a", 
                  "CD154" = "CD154", "IL2" = "IL-2", "IL17a" = "IL-17a", 
                  "IL4/5/13" = "IL-4/5/13", "IFNg" = "IFN-\u03B3", 
                  "CCR7" = "CCR7", "CD45RA" = "CD45RA",
                  "CD38" = "CD38", "HLADR" = "HLA-DR")
plot_scaled_mfi_v2 <- function(currentColumn) {
  myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
  sc_mfi <- scale_colour_gradientn(colours = myPalette(100), oob=squish, limits=c(-2.5, 2.5))
  gg_themes +
    geom_point_rast(aes(colour=cd4_compass_subsets_w_umap[,paste0(currentColumn, ".scaled")]),
               shape=20, alpha=0.8, size=pointSize) + # !!as.name(paste0(currentColumn, ".scaled")) doesn't work
    sc_mfi +
    ggtitle(mfi_col_text_4plot[[currentColumn]])
}

unstratified_mfi_plot_cols <- c("TNFa", "CD107a", 
                  "CD154", "IL2", "IL17a", 
                  "IL4/5/13", "IFNg", 
                  "CCR7", "CD45RA",
                  "CD38", "HLADR")
mfi_plots_unstrat <- lapply(unstratified_mfi_plot_cols, plot_scaled_mfi_v2)
names(mfi_plots_unstrat) <- unstratified_mfi_plot_cols
# Memory UMAP plot colored by gate membership (TEMRA, TEM, TCM, Naive)
getMemoryCategory <- function(CD45RA, CCR7) {
  if(CD45RA & CCR7) {
    "Naive"
  } else if(CD45RA & !CCR7) {
    "TEMRA"
  } else if(!CD45RA & CCR7) {
    "TCM"
  } else if(!CD45RA & !CCR7) {
    "TEM"
  } else {
    "Uncategorized"
  }
}

cd4_compass_subsets_w_umap$MemoryCategory <- pmap_chr(.l = list(cd4_compass_subsets_w_umap$`/Time/LD-3+/1419-3+/S/Lymph/4+/CD45RA+`,
                                                               cd4_compass_subsets_w_umap$`/Time/LD-3+/1419-3+/S/Lymph/4+/CCR7+`),
                                                     .f = getMemoryCategory)

myPalette(4)
## [1] "#5E4FA2" "#BEE4A0" "#FDBE6F" "#9E0142"
scales::show_col(myPalette(4))

memColorScheme <- c("TEM" = "#9E0142", "TEMRA" = "#BEE4A0", "TCM" = "#5E4FA2", "Naive" = "#FDBE6F")
mem_plot_colored_by_gate <- gg_themes +
  geom_point_rast(aes(colour=cd4_compass_subsets_w_umap[,"MemoryCategory"]),
               shape=20, alpha=0.8, size=pointSize) +
  ggtitle("Memory") +
  scale_color_manual(values = memColorScheme)
# + theme(legend.position = "right") + guides(colour = guide_legend(override.aes = list(size=9)))

mem_plot_colored_by_gate

# Activation plots: color by gate membership, one plot per marker
boolColorScheme <- c("0" = "#E2E2E2", "1" = "#363636")
hladr_bool_plot <- gg_themes +
  geom_point_rast(aes(colour=factor(cd4_compass_subsets_w_umap[,"/Time/LD-3+/1419-3+/S/Lymph/HLADR+"], levels = c(1, 0))),
               shape=20, alpha=0.8, size=pointSize) +
  scale_color_manual(values = boolColorScheme, labels=c("1"="Expressed", "0"="Not Expressed")) +
  ggtitle("HLA-DR")
cd38_bool_plot <- gg_themes +
  geom_point_rast(aes(colour=factor(cd4_compass_subsets_w_umap[,"/Time/LD-3+/1419-3+/S/Lymph/CD38+"], levels = c(1, 0))),
               shape=20, alpha=0.8, size=pointSize) +
  scale_color_manual(values = boolColorScheme, labels=c("1"="Expressed", "0"="Not Expressed")) +
  ggtitle("CD38")
# factor_colors = hsv((seq(0, 1, length.out = 4 + 1)[-1] +
#                          0.2)%%1, 0.7, 0.95)
# 
# stim_colors <- c("Spike 1" = "#fc68be", "Spike 2" = "#d94e09", "NCAP" = "#6B49F2", "VEMP" = "#D0F249", "DMSO" = "#91570a") # heatmap colors
# stim_colors <- c("Spike 1" = "#49F2BF", "Spike 2" = "#F2497C", "NCAP" = "#6B49F2", "VEMP" = "#D0F249", "DMSO" = "#91570a")
stim_abbreviations = c("Spike 1" = "S1", "Spike 2" = "S2", "NCAP" = "NCAP", "VEMP" = "VEMP", "DMSO" = "DMSO")

plot_stim_contour_plot <- function(currentStim) {
  gg_themes +
    geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
    geom_density_2d(data=cd4_compass_subsets_w_umap %>% dplyr::filter(STIM == currentStim),
                    colour="#363636", bins=12) + # stim_colors[[currentStim]]
    ggtitle(stim_abbreviations[[currentStim]])
}

stims <- c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")
stim_plots <- lapply(stims, plot_stim_contour_plot)
names(stim_plots) <- stims
# Unstratified PolyF plot
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
sc_polyf <- scale_colour_gradientn(colours = myPalette(100), limits=c(1, 6))
polyf_plot_unstrat <- gg_themes +
  geom_point_rast(aes(colour=cd4_compass_subsets_w_umap[,"cytokine_degree"]), shape=20, alpha=0.8, size=pointSize) +
  sc_polyf +
  ggtitle("PolyF")
# Extract legends
mfi_legend <- as_ggplot(get_legend(mfi_plots_unstrat$TNFa +
                                     theme(legend.position="bottom",
                                             legend.title.align = 0.5) +
                                     labs(colour="Scaled Expression") +
                                     guides(colour=guide_colorbar(title.position = "top"))))
polyf_legend <- as_ggplot(get_legend(polyf_plot_unstrat +
                                       theme(legend.position="bottom",
                                             legend.title.align = 0.5) +
                                       labs(colour="PolyF") +
                                       guides(colour=guide_colorbar(title.position = "top"))))
mem_legend <- as_ggplot(get_legend(mem_plot_colored_by_gate +
                                     theme(legend.position = "right",
                                           legend.title.align = 0.5) +
                                     guides(colour = guide_legend(override.aes = list(size=9))) +
                                     labs(colour="Memory")))
bool_legend <- as_ggplot(get_legend(hladr_bool_plot +
                                     theme(legend.position = "right",
                                           legend.title.align = 0.5) +
                                     guides(colour = guide_legend(override.aes = list(size=9))) +
                                     labs(colour="Gate Membership")))


mfi_legend_vertical <- as_ggplot(get_legend(mfi_plots_unstrat$TNFa +
                                     theme(legend.position="right",
                                           legend.title.align = 0.5,
                                           legend.title = element_text(angle = 90)) +
                                     labs(colour="Scaled Expression") +
                                     guides(colour=guide_colorbar(title.position = "left"))))
polyf_legend_vertical <- as_ggplot(get_legend(polyf_plot_unstrat +
                                       theme(legend.position="right",
                                             legend.title.align = 0.5,
                                             legend.title = element_text(angle = 90)) +
                                       labs(colour="PolyF") +
                                       guides(colour=guide_colorbar(title.position = "left"))))

Put it all together

print(
  # Top row
  (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
  stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
  polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
  # Second row
  (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
  stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
  mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
  # Third row
  (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
  stim_plots$DMSO | plot_spacer() | plot_spacer() |
  mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
)

if(save_output) {
  png(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP.png",
                              date)),
      width=18, height=9, units="in", res=300)
    print(
      # Top row
      (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
      stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
      polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
      # Second row
      (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
      stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
      mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
      # Third row
      (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
      stim_plots$DMSO | plot_spacer() | plot_spacer() |
      mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
    )
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_MFI_Legend.png",
                            date)),
    width=3, height=1.5, units="in", res=300)
  print(mfi_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_PolyF_Legend.png",
                            date)),
    width=3, height=1.5, units="in", res=300)
  print(polyf_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_Memory_Legend.png",
                            date)),
    width=2, height=3, units="in", res=300)
  print(mem_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_GateMembership_Legend.png",
                          date)),
  width=2, height=3, units="in", res=300)
  print(bool_legend)
  dev.off()
  
  # PDFs
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_rast.pdf",
                              date)),
      width=18, height=9, onefile = TRUE, bg = "transparent", family = "Arial")
  print(
      # Top row
      (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
      stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
      polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
      # Second row
      (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
      stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
      mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
      # Third row
      (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
      stim_plots$DMSO | plot_spacer() | plot_spacer() |
      mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
    )
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_MFI_Legend.pdf",
                            date)),
    width=3, height=1.5, onefile = TRUE, bg = "transparent", family = "Arial")
  print(mfi_legend)
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_PolyF_Legend.pdf",
                            date)),
    width=3, height=1.5, onefile = TRUE, bg = "transparent", family = "Arial")
  print(polyf_legend)
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_Memory_Legend.pdf",
                            date)),
    width=2, height=3, onefile = TRUE, bg = "transparent", family = "Arial")
  print(mem_legend)
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_GateMembership_Legend.pdf",
                          date)),
  width=2, height=3, onefile = TRUE, bg = "transparent", family = "Arial")
  print(bool_legend)
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_MFI_Legend_Vertical.pdf",
                            date)),
    width=1.5, height=3, onefile = TRUE, bg = "transparent", family = "Arial")
  print(mfi_legend_vertical)
  dev.off()
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD4_COMPASS_Subsets_UMAP_PolyF_Legend_Vertical.pdf",
                            date)),
    width=1.5, height=3, onefile = TRUE, bg = "transparent", family = "Arial")
  print(polyf_legend_vertical)
  dev.off()
}

CD8

# Drop certain wells for just the CD8 runs due to low CD8 count:
gs3 <- subset(gs2, !(STIM %in% c("Spike 2", "NCAP") & `SAMPLE ID` %in% c("BWT20", "15548")) &
         !(STIM == "Spike 2" & `SAMPLE ID` == "15530"))

Extract mfi data

cd8_gates_for_dimred <- c(
  "/Time/LD-3+/1419-3+/S/Lymph/8+/107a", "/Time/LD-3+/1419-3+/S/Lymph/8+/154", 
  "/Time/LD-3+/1419-3+/S/Lymph/8+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL2", 
  "/Time/LD-3+/1419-3+/S/Lymph/8+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL4513", 
  "/Time/LD-3+/1419-3+/S/Lymph/8+/TNF",
  "/Time/LD-3+/1419-3+/S/Lymph/8+/CCR7+", "/Time/LD-3+/1419-3+/S/Lymph/8+/CD45RA+",
  "/Time/LD-3+/1419-3+/S/Lymph/CD38+", "/Time/LD-3+/1419-3+/S/Lymph/HLADR+")
cd8_cytokine_gates <- c("/Time/LD-3+/1419-3+/S/Lymph/8+/107a", "/Time/LD-3+/1419-3+/S/Lymph/8+/154", 
                        "/Time/LD-3+/1419-3+/S/Lymph/8+/IFNG", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL2", 
                        "/Time/LD-3+/1419-3+/S/Lymph/8+/IL17", "/Time/LD-3+/1419-3+/S/Lymph/8+/IL4513", 
                        "/Time/LD-3+/1419-3+/S/Lymph/8+/TNF")
cd8_compass_subsets_parentGate <- "/Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets"

pop_counts <- pData(gs3) %>% 
  left_join(gs_pop_get_count_fast(gs3, subpopulations = cd8_compass_subsets_parentGate),
            by = c("rowname" = "name")) %>% 
  dplyr::rename(CD8_COMPASS_Subsets = Count) %>% 
  dplyr::select(rowname, Batch, "SAMPLE ID", STIM, Cohort, CD8_COMPASS_Subsets) %>% 
  dplyr::filter(!(Cohort %in% c(NA, "Healthy control", "Healthy control 2017-2018")) & STIM != "SEB")

Keep in mind that there is lopsided patient and group representation simply due to not sampling:

cd8_compass_subsets_sampleSizes_4plot <- pop_counts %>%
  mutate(Cohort = factor(Cohort,
                         levels = c("Non-hospitalized", "Hospitalized"),
                         labels = c("Conv\nNon-Hosp", "Conv\nHosp")))
ggplot(cd8_compass_subsets_sampleSizes_4plot,
       aes(factor(Cohort), CD8_COMPASS_Subsets)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.15, height = 0) +
  theme_bw(base_size=20) +
  labs(title="ICS CD8 UMAP patient representation",
       y="CD8+ COMPASS Subset+ Events\n for Dimensionality Reduction\n(not sampled)") +
  facet_grid(Batch ~ STIM) +
  theme(axis.title.x = element_blank())

# Extract data for dimensionality reduction (not actually sampling)
call_sampleGatingHierarchy_for_cd8 <- function(currentSampleName) {
  # print(sprintf("Sampling data from %s", currentSampleName))
  sampleGatingHierarchy(gs3[[currentSampleName]], cd8_compass_subsets_parentGate, n = NULL, otherGates = cd8_gates_for_dimred)
}
cd8_compass_subsets_data <- map_dfr(pop_counts$rowname, call_sampleGatingHierarchy_for_cd8)
dim(cd8_compass_subsets_data)
## [1] 19904    54
knitr::kable(head(cd8_compass_subsets_data))
Days symptom onset to visit 1 Sex WELL ID Pair ID Hispanic? STIM name Cell count Race Cohort Age Race_v2 rowname $DATE Batch Sample ID SAMPLE ID PLATE NAME Collection date EXPERIMENT NAME filename /Time/LD-3+/1419-3+/S/Lymph/8+/CD8_COMPASS_Subsets /Time/LD-3+/1419-3+/S/Lymph/8+/107a /Time/LD-3+/1419-3+/S/Lymph/8+/154 /Time/LD-3+/1419-3+/S/Lymph/8+/IFNG /Time/LD-3+/1419-3+/S/Lymph/8+/IL2 /Time/LD-3+/1419-3+/S/Lymph/8+/IL17 /Time/LD-3+/1419-3+/S/Lymph/8+/IL4513 /Time/LD-3+/1419-3+/S/Lymph/8+/TNF /Time/LD-3+/1419-3+/S/Lymph/8+/CCR7+ /Time/LD-3+/1419-3+/S/Lymph/8+/CD45RA+ /Time/LD-3+/1419-3+/S/Lymph/CD38+ /Time/LD-3+/1419-3+/S/Lymph/HLADR+ Time FSC-A FSC-H SSC-A SSC-H CD8b TNFa CD107a CD154 CD3 ECD IL2 CD4 IL17a IL4/5/13 CD14/CD19 CCR7 CD38 L/D IFNg CD45RA HLADR
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 1 0 1 0 0 0 0 0 1 0 0 15.665 93581.08 76447 20852.16 20555 1791.021 774.6697 2354.91577 763.7285 2899.465 757.1143 1151.5387 463.08521 462.8966 212.82829 1400.210 964.2203 1244.181 1641.4823 2330.406 304.0638
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 1 0 1 1 0 0 16.102 137163.28 119576 16784.91 16232 1802.571 391.5926 377.20135 1128.8015 2867.831 802.3138 1174.3997 1619.48035 2986.4995 1120.99817 2181.621 1053.1833 1177.822 736.3119 2701.324 772.7253
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 1 0 0 0 0 0 0 1 1 0 0 16.680 89083.40 73369 23013.24 22141 1625.985 424.0639 1784.99561 316.6937 3008.369 865.6202 419.3799 444.15247 975.1948 11.76863 2401.309 1482.0504 1132.399 507.0199 2737.020 775.5408
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 1 0 1 1 0 0 18.756 111505.68 96285 18827.67 18237 2493.273 597.3187 575.55695 986.9149 3061.338 953.5745 987.6654 77.79742 1512.9031 1129.35010 2317.092 1527.4861 1221.771 934.6315 2658.387 1019.6381
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 0 0 0 1 0 1 1 0 0 18.822 136115.23 121886 18064.68 18094 1878.481 624.7535 504.89352 1381.5298 3008.781 451.4085 692.4930 1031.82068 1667.1205 865.71771 2200.889 1423.8490 1090.599 526.6213 2925.936 717.5511
61 M H12 11 N NCAP 112590.fcs 1.6x10^7 White Non-hospitalized 33 White 112590.fcs_366900 28-MAY-2020 1 90C.1.A 90C P2 2020-05-01 20200528_COVID_ICS-B1 90C_H12_H12_096.fcs 1 0 0 1 0 0 0 0 1 1 0 0 20.038 105345.12 89530 28575.15 26810 2319.511 174.2149 -19.19797 690.6313 3158.891 431.6187 1156.5277 537.53528 703.9210 908.09760 2628.847 668.0073 1490.596 1904.6018 2970.056 911.3188

Is there maybe one cytokine that is dominating the entire sample?
If one cytokine has very high background expression (and a generous gate), it could be gated positive in a lot of events.
The high number of events expressing this cytokine could lead to it dominating the data, so that most sampled events are positive for this noisy cytokine. It would drown out real signal from other cytokines.

cytokine_dominance <- cd8_compass_subsets_data %>% 
  group_by(Batch) %>% 
  summarise_at(cd8_cytokine_gates, sum) %>%
  t() %>%
  as.data.frame() %>%
  set_names(c("B1", "B2", "B3")) %>%
  rownames_to_column("Cytokine_Gate") %>%
  dplyr::filter(Cytokine_Gate != "Batch")
knitr::kable(cytokine_dominance)
Cytokine_Gate B1 B2 B3
/Time/LD-3+/1419-3+/S/Lymph/8+/107a 2008 5318 3449
/Time/LD-3+/1419-3+/S/Lymph/8+/154 0 0 0
/Time/LD-3+/1419-3+/S/Lymph/8+/IFNG 666 939 609
/Time/LD-3+/1419-3+/S/Lymph/8+/IL2 56 119 59
/Time/LD-3+/1419-3+/S/Lymph/8+/IL17 0 0 0
/Time/LD-3+/1419-3+/S/Lymph/8+/IL4513 1511 3568 2405
/Time/LD-3+/1419-3+/S/Lymph/8+/TNF 30 341 140
cytokine_dominance %>%
  pivot_longer(cols = starts_with("B"), names_to = "Batch", values_to = "Events_in_Gate") %>% 
  mutate(Cytokine = sub(".*4\\+\\/(.*)", "\\1", Cytokine_Gate)) %>% 
  ggplot(aes(Cytokine, Events_in_Gate, fill = Cytokine)) +
  theme_bw(base_size=18) +
  geom_bar(stat="identity") +
  facet_grid(. ~ Batch) +
  labs(title = "CD8 Run Cytokine Dominance by Batch")

Perform Dimensionality Reduction

cols_4_dimred <- c("CD3 ECD", "CD8b", "CD4",
                  "TNFa", "CD107a", 
                  "CD154", "IL2", "IL17a", 
                  "IL4/5/13", "IFNg", 
                  "CCR7", "CD45RA",
                  "CD38", "HLADR")
cd4.scaled_dimred_input <- cd8_compass_subsets_data %>% 
  dplyr::select(Batch, all_of(cols_4_dimred)) %>% 
  group_by(Batch) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(data = lapply(data, function(df) {as.data.frame(scale(as.matrix(df)))})) %>% 
  unnest(cols = c(data)) %>% 
  rename_at(vars(all_of(cols_4_dimred)),function(x) paste0(x,".scaled")) %>% 
  dplyr::select(-Batch)

cd8_compass_subsets_data <- cbind(cd8_compass_subsets_data, cd4.scaled_dimred_input)

# UMAP can take a long time, so there is a rerun_dimred switch
if(rerun_dimred) {
  print("Running UMAP")
  set.seed(date)
  print(Sys.time())
  cd8_compass_subsets_dimred_out <- cd8_compass_subsets_data %>% 
    # Run CD3, co-receptor, cytokine, memory, and activation markers through UMAP
    dplyr::select(all_of(paste0(cols_4_dimred, ".scaled"))) %>% 
    uwot::umap(spread = 9, min_dist = 0.02, n_threads = 7)
  print(Sys.time())
  cd8_compass_subsets_w_umap <- cbind(as.data.frame(cd8_compass_subsets_dimred_out) %>% 
    dplyr::rename(x.umap = V1, y.umap = V2),
    cd8_compass_subsets_data)
  if(save_output) {
    saveRDS(cd8_compass_subsets_w_umap, here::here(sprintf("out/UMAP/%s_ICS_CD8_COMPASS_Subsets_UMAP_Unsampled.rds", date)))
  }
} else {
  # Assuming UMAP results are already saved
  print("Loading saved UMAP run")
  cd8_compass_subsets_w_umap <- readRDS(here::here(sprintf("out/UMAP/%s_ICS_CD8_COMPASS_Subsets_UMAP_Unsampled.rds", date)))
}
## [1] "Loading saved UMAP run"

Plot UMAP results

Shuffle data frame rows so e.g. Batch 3 doesn’t dominate foreground

set.seed(date)
cd8_compass_subsets_w_umap <- cd8_compass_subsets_w_umap[sample(nrow(cd8_compass_subsets_w_umap), nrow(cd8_compass_subsets_w_umap)),]
# Arial font setup. Downloaded afms from https://github.com/microsoft/microsoft-r-open/tree/ec3fd89e5fb5794bd8149905c134ad801bb61800
Arial <- Type1Font(family = "Arial",
                   metrics = c(here::here("data/Arial_afm/ArialMT.afm"),
                               here::here("data/Arial_afm/ArialMT-Bold.afm"), 
                               here::here("data/Arial_afm/ArialMT-Italic.afm"),
                               here::here("data/Arial_afm/ArialMT-BoldItalic.afm")))
pdfFonts(Arial = Arial)

boolColorScheme <- c("FALSE" = "#E2E2E2", "TRUE" = "#023FA5")

cd8_compass_subsets_w_umap <- cd8_compass_subsets_w_umap %>% 
  mutate(cytokine_degree = rowSums(dplyr::select(., all_of(cd8_cytokine_gates))))
cd8_compass_subsets_w_umap <- cd8_compass_subsets_w_umap %>% 
  mutate(Cohort = factor(Cohort, levels = c("Non-hospitalized", "Hospitalized")),
         STIM = factor(STIM, levels = c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")))

stim_labs <- c("DMSO", "S1", "S2", "NCAP", "VEMP") 
names(stim_labs) <- c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")
cohort_labs <- c("Conv\nNon-Hosp", "Conv\nHosp")
names(cohort_labs) <- c("Non-hospitalized", "Hospitalized")

base_dimred_plot <- function(currentColumn, pointSize = 0.02, colorScheme = NA) {
  p <- ggplot(cd8_compass_subsets_w_umap, aes(x=x.umap, y=y.umap,
                                       colour=if(currentColumn %in% c("Batch", "SAMPLE ID")) {
                                         factor(!!as.name(currentColumn))
                                         } else {
                                           as.logical(!!as.name(currentColumn))
                                           })) +
    geom_point_rast(shape=20, alpha=0.8, size=pointSize) +
    facet_grid(Cohort ~ STIM, switch="y",
               labeller = labeller(Cohort = cohort_labs, STIM = stim_labs)) +
    theme_bw() +
    theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
          axis.text=element_blank(),
          axis.line=element_blank(),
          axis.ticks=element_blank(),
          axis.title=element_blank(),
          strip.text=element_text(face="bold", size=22),
          panel.grid.major = element_blank(),
          legend.title=element_text(face="bold", size=14),
          strip.text.x = element_text(margin = margin(0.15,0,0.15,0, "cm")))
  if(!anyNA(colorScheme)) {
    p <- p + scale_color_manual(values = colorScheme)
  }
  if(currentColumn %in% c("Batch", "SAMPLE ID")) {
    p <- p + labs(color = currentColumn) +
      guides(colour = guide_legend(override.aes = list(size=10))) +
      theme(legend.title = element_text(size=15),
            legend.text = element_text(size=15),
            legend.position = "bottom") +
      scale_colour_manual(values=as.character(iwanthue(length(unique(cd8_compass_subsets_w_umap[,currentColumn])))))
  } else {
    p <- p + theme(legend.position = "none")
  }
  p
}

Look for Batch and Patient sub-clustering

base_dimred_plot("Batch")

base_dimred_plot("SAMPLE ID")

Cytokine, Memory, and Activation expression localization

Now visualize the cytokine, memory, and activation expression localization

for(cg in cd8_gates_for_dimred) {
  print(base_dimred_plot(cg, colorScheme = boolColorScheme) +
    labs(title = sprintf("CD8+ COMPASS Subset+ UMAP\nColored by %s", sub(".*\\/([^(\\/)]+)", "\\1", cg))))
}

Color by degree

table(cd8_compass_subsets_w_umap$cytokine_degree)
## 
##     1     2     3 
## 18861   772   271
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
sc <- scale_colour_gradientn(colours = myPalette(100), limits=c(1, 6))

currentColumn <- "cytokine_degree"
pointSize <- 0.02
ggplot(cd8_compass_subsets_w_umap, aes(x=x.umap, y=y.umap, colour=!!as.name(currentColumn))) +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize) +
  facet_grid(Cohort ~ STIM, switch="y",
             labeller = labeller(Cohort = cohort_labs, STIM = stim_labs)) +
  labs(colour="Cytokine\nDegree",
       title="CD8+ COMPASS Subset+ UMAP\nColored by Cytokine Polyfunctionality") +
  theme_bw() +
  theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        strip.text=element_text(face="bold", size=22),
        panel.grid.major = element_blank(),
        legend.title=element_text(face="bold", size=15),
        legend.text = element_text(size=13),
        strip.text.x = element_text(margin = margin(0.15,0,0.15,0, "cm")),
        legend.position = "bottom") +
  sc

Unstratified UMAP plots

pointSize <- 0.02
# Theme settings
# Expand axis limits so contour lines don't get cut off near borders
gg_themes <- ggplot(cd8_compass_subsets_w_umap, aes(x=x.umap, y=y.umap)) +
  scale_x_continuous(limits=range(cd8_compass_subsets_w_umap$x.umap)*1.16) +
  scale_y_continuous(limits=range(cd8_compass_subsets_w_umap$y.umap)*1.11) +
  theme_bw() +
  theme(plot.title=element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = margin(1, 0, 0, 0, unit = "pt"))
# Hospitalization status contour plots
hosp_contour_plot <- gg_themes +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
  geom_density_2d(data=cd8_compass_subsets_w_umap %>% dplyr::filter(Cohort == "Hospitalized"),
                  colour="#363636", bins=12) + # "red"
  ggtitle("Hosp")

nonhosp_contour_plot <- gg_themes +
  geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
  geom_density_2d(data=cd8_compass_subsets_w_umap %>% dplyr::filter(Cohort == "Non-hospitalized"),
                  colour="#363636", bins=12) + # "blue"
  ggtitle("Not-Hosp")
# Scaled mfi expression plots for memory, activation, and cytokine markers 
mfi_col_text_4plot <- c("TNFa" = "TNF", "CD107a" = "CD107a", 
                  "CD154" = "CD154", "IL2" = "IL-2", "IL17a" = "IL-17a", 
                  "IL4/5/13" = "IL-4/5/13", "IFNg" = "IFN-\u03B3", 
                  "CCR7" = "CCR7", "CD45RA" = "CD45RA",
                  "CD38" = "CD38", "HLADR" = "HLA-DR")
plot_scaled_mfi_v2 <- function(currentColumn) {
  myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
  sc_mfi <- scale_colour_gradientn(colours = myPalette(100), oob=squish, limits=c(-2.5, 2.5))
  gg_themes +
    geom_point_rast(aes(colour=cd8_compass_subsets_w_umap[,paste0(currentColumn, ".scaled")]),
               shape=20, alpha=0.8, size=pointSize) + # !!as.name(paste0(currentColumn, ".scaled")) doesn't work
    sc_mfi +
    ggtitle(mfi_col_text_4plot[[currentColumn]])
}

unstratified_mfi_plot_cols <- c("TNFa", "CD107a", 
                  "CD154", "IL2", "IL17a", 
                  "IL4/5/13", "IFNg", 
                  "CCR7", "CD45RA",
                  "CD38", "HLADR")
mfi_plots_unstrat <- lapply(unstratified_mfi_plot_cols, plot_scaled_mfi_v2)
names(mfi_plots_unstrat) <- unstratified_mfi_plot_cols
# Memory UMAP plot colored by gate membership (TEMRA, TEM, TCM, Naive)
getMemoryCategory <- function(CD45RA, CCR7) {
  if(CD45RA & CCR7) {
    "Naive"
  } else if(CD45RA & !CCR7) {
    "TEMRA"
  } else if(!CD45RA & CCR7) {
    "TCM"
  } else if(!CD45RA & !CCR7) {
    "TEM"
  } else {
    "Uncategorized"
  }
}

cd8_compass_subsets_w_umap$MemoryCategory <- pmap_chr(.l = list(cd8_compass_subsets_w_umap$`/Time/LD-3+/1419-3+/S/Lymph/8+/CD45RA+`,
                                                               cd8_compass_subsets_w_umap$`/Time/LD-3+/1419-3+/S/Lymph/8+/CCR7+`),
                                                     .f = getMemoryCategory)

memColorScheme <- c("TEM" = "#9E0142", "TEMRA" = "#BEE4A0", "TCM" = "#5E4FA2", "Naive" = "#FDBE6F")
mem_plot_colored_by_gate <- gg_themes +
  geom_point_rast(aes(colour=cd8_compass_subsets_w_umap[,"MemoryCategory"]),
               shape=20, alpha=0.8, size=pointSize) +
  ggtitle("Memory") +
  scale_color_manual(values = memColorScheme)
# + theme(legend.position = "right") + guides(colour = guide_legend(override.aes = list(size=9)))
# Activation plots: color by gate membership, one plot per marker
boolColorScheme <- c("0" = "#E2E2E2", "1" = "#363636")
hladr_bool_plot <- gg_themes +
  geom_point_rast(aes(colour=factor(cd8_compass_subsets_w_umap[,"/Time/LD-3+/1419-3+/S/Lymph/HLADR+"], levels = c(1, 0))),
               shape=20, alpha=0.8, size=pointSize) +
  scale_color_manual(values = boolColorScheme, labels=c("1"="Expressed", "0"="Not Expressed")) +
  ggtitle("HLA-DR")
cd38_bool_plot <- gg_themes +
  geom_point_rast(aes(colour=factor(cd8_compass_subsets_w_umap[,"/Time/LD-3+/1419-3+/S/Lymph/CD38+"], levels = c(1, 0))),
               shape=20, alpha=0.8, size=pointSize) +
  scale_color_manual(values = boolColorScheme, labels=c("1"="Expressed", "0"="Not Expressed")) +
  ggtitle("CD38")
# factor_colors = hsv((seq(0, 1, length.out = 4 + 1)[-1] +
#                          0.2)%%1, 0.7, 0.95)
# 
# stim_colors <- c("Spike 1" = "#fc68be", "Spike 2" = "#d94e09", "NCAP" = "#6B49F2", "VEMP" = "#D0F249", "DMSO" = "#91570a") # heatmap colors
# stim_colors <- c("Spike 1" = "#49F2BF", "Spike 2" = "#F2497C", "NCAP" = "#6B49F2", "VEMP" = "#D0F249", "DMSO" = "#91570a")
stim_abbreviations = c("Spike 1" = "S1", "Spike 2" = "S2", "NCAP" = "NCAP", "VEMP" = "VEMP", "DMSO" = "DMSO")

plot_stim_contour_plot <- function(currentStim) {
  gg_themes +
    geom_point_rast(shape=20, alpha=0.8, size=pointSize, colour="#E2E2E2") +
    geom_density_2d(data=cd8_compass_subsets_w_umap %>% dplyr::filter(STIM == currentStim),
                    colour="#363636", bins=12) + # stim_colors[[currentStim]]
    ggtitle(stim_abbreviations[[currentStim]])
}

stims <- c("DMSO", "Spike 1", "Spike 2", "NCAP", "VEMP")
stim_plots <- lapply(stims, plot_stim_contour_plot)
names(stim_plots) <- stims
# Unstratified PolyF plot
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
sc_polyf <- scale_colour_gradientn(colours = myPalette(100), limits=c(1, 6))
polyf_plot_unstrat <- gg_themes +
  geom_point_rast(aes(colour=cd8_compass_subsets_w_umap[,"cytokine_degree"]), shape=20, alpha=0.8, size=pointSize) +
  sc_polyf +
  ggtitle("PolyF")
# Extract legends
mfi_legend <- as_ggplot(get_legend(mfi_plots_unstrat$TNFa +
                                     theme(legend.position="bottom",
                                             legend.title.align = 0.5) +
                                     labs(colour="Scaled Expression") +
                                     guides(colour=guide_colorbar(title.position = "top"))))
polyf_legend <- as_ggplot(get_legend(polyf_plot_unstrat +
                                       theme(legend.position="bottom",
                                             legend.title.align = 0.5) +
                                       labs(colour="PolyF") +
                                       guides(colour=guide_colorbar(title.position = "top"))))
mem_legend <- as_ggplot(get_legend(mem_plot_colored_by_gate +
                                     theme(legend.position = "right",
                                           legend.title.align = 0.5) +
                                     guides(colour = guide_legend(override.aes = list(size=9))) +
                                     labs(colour="Memory")))
bool_legend <- as_ggplot(get_legend(hladr_bool_plot +
                                     theme(legend.position = "right",
                                           legend.title.align = 0.5) +
                                     guides(colour = guide_legend(override.aes = list(size=9))) +
                                     labs(colour="Gate Membership")))

Put it all together

print(
  # Top row
  (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
  stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
  polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
  # Second row
  (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
  stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
  mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
  # Third row
  (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
  stim_plots$DMSO | plot_spacer() | plot_spacer() |
  mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
)

if(save_output) {
  png(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP.png",
                              date)),
      width=18, height=9, units="in", res=300)
    print(
      # Top row
      (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
      stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
      polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
      # Second row
      (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
      stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
      mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
      # Third row
      (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
      stim_plots$DMSO | plot_spacer() | plot_spacer() |
      mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
    )
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP_MFI_Legend.png",
                            date)),
    width=3, height=1.5, units="in", res=300)
  print(mfi_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP_PolyF_Legend.png",
                            date)),
    width=3, height=1.5, units="in", res=300)
  print(polyf_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP_Memory_Legend.png",
                            date)),
    width=2, height=3, units="in", res=300)
  print(mem_legend)
  dev.off()
  
  png(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP_GateMembership_Legend.png",
                          date)),
  width=2, height=3, units="in", res=300)
  print(bool_legend)
  dev.off()
  
  # PDFs
  
  cairo_pdf(file=here::here(sprintf("out/UMAP/%s_CD8_COMPASS_Subsets_UMAP_rast.pdf",
                            date)),
    width=18, height=9, onefile = TRUE, bg = "transparent", family = "Arial")
    print(
      # Top row
      (hosp_contour_plot | nonhosp_contour_plot | plot_spacer() |
      stim_plots$`Spike 1` | stim_plots$`Spike 2` | plot_spacer() |
      polyf_plot_unstrat | plot_spacer() | mfi_plots_unstrat$IFNg) /
      # Second row
      (mem_plot_colored_by_gate | plot_spacer() | plot_spacer() |
      stim_plots$NCAP | stim_plots$VEMP | plot_spacer() |
      mfi_plots_unstrat$TNFa | mfi_plots_unstrat$IL2 | mfi_plots_unstrat$`IL4/5/13`) /
      # Third row
      (hladr_bool_plot | cd38_bool_plot | plot_spacer() |
      stim_plots$DMSO | plot_spacer() | plot_spacer() |
      mfi_plots_unstrat$IL17a | mfi_plots_unstrat$CD107a | mfi_plots_unstrat$CD154 )
    )
  dev.off()
}